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We investigate geometric resonances in Bose-Einstein condensates by solving the underlying time- 
dependent Gross-Pitaevskii (GP) equation for systems with two-and three-body interactions in an 
axially-symmetric harmonic trap. To this end, we use a recently developed analytical method 
[Phys. Rev. A 84, 013618 (2011)], based on a perturbative expansion and a Poincare-Lindstedt 
analysis of a Gaussian variational approach, as well as a detailed numerical study of a set of ordi- 
nary differential equations for variational parameters. By changing the anisotropy of the confining 
potential, we numerically observe and analytically describe strong nonlinear effects: shifts in the 
frequencies and mode coupling of collective modes, as well as resonances. Furthermore, we discuss in 
detail the stability of a Bose-Einstein condensate in the presence of an attractive two-body interac- 
tion and a repulsive three-body interaction. We show that a small repulsive three-body interaction 
is able to extend the stability region of the condensate as it increases the critical number of atoms 
in the trap. 

PACS numbers: 03.75.Kk, 03.75.Nt, 67.85. De 



I. INTRODUCTION 



The experimental discovery of Bose-Einstein condensa- 
tion in dilute atomic vapors [11-01 has instigated a broad 
interest in ultracold atoms and molecules, and paved the 
way for extensive studies of a wide range of experimen- 
tal and theoretical topics. In particular, many exper- 
iments have focused on investigating collective excita- 
tions of harmonically trapped Bose-Einstein condensates 
(BECs), as they can be measured very accurately and 
thus provide a reliable method for extracting the respec- 
tive system parameters of such ultracold systems 0]. 

The excellent agreement between the measured fre- 
quencies [sl-fTol and theoretical predictions [TTl - [l7j is one 
of the first important achievements in the investiga- 
tion of such systems. At the mean- field level, they can 
be successfully described by the time-dependent Gross- 
Piteavskii equation for the macroscopic wave function of 
a BEC at zero temperature (TB-HH • Due to the inherent 
nonlinearity in this equation of motion, a wide variety of 
interesting phenomena are observed in collective excita- 
tions of BECs, including frequency shifts [22l . [23| , mode 
couphn g |2^ . [2^ |25| , damping 0, Ull , nonlinear interfer- 
ometry 1^ , as well as collapse and revival of oscillations 
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The collective oscillation modes can be induced in a 
BEC by modulating the external potential trap 
[22l . IsOmSt . Alternatively, this can also be allowed by a 
modulation of the s-wave scattering length [23;, ii6,- 38] or 
three-body interactions 0, |3^. In the Thomas- Fermi 
(TF) limit, Stringari [ll| has analytically calculated fre- 
quencies of the collective modes, while Edwards et al. 
[12| extended these calculations via a numerical solution 
of the GP equation. 

BEC systems are highly nonlinear and, therefore, a 
resonant coupling between collective modes is expected. 
It was experimentally observed [13, E^l that, when the 
parity quadrupole mode is excited by changing the trap 
anisotropy parameter above a certain value, it is possible 
to achieve an energy transfer between modes at a rate 
which is comparable to the collective mode frequency. In 
Ref. 22| , the frequency shift of collective modes was ana- 
lytically studied due to the change of the trap anisotropy 
in a generic axially-symmetric geometry using the hydro- 
dynamic equations in the Thomas-Fermi approximation 
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In this way it was shown that the collective modes 
exhibit a resonant behavior for specific values of the trap 
anisotropy, which are called geometric resonances, so that 
strong nonlinear effects can be observed even for oscilla- 
tions of relatively small amplitude. In Ref. [4l|, the au- 
thors studied the dynamics of large amplitude collective 
mode oscillations, frequency shifts and mode couplings 
for the case of a superfluid Fermi gas in the transition 
from a BCS superfluid to a BEC based on a superfluid hy- 
drodynamic approach. The coupling rates between vari- 
ous modes were calculated analytically in Ref. [11], and 
it was shown that the coupling can be well described by 
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a simple Hamiltonian, enabling quantitative studies of 
the squeezing effects which are related to the harmonic- 
generation processes. The excitations of quadrupole and 
scissors modes in two-component BECs with a particular 
emphasis on the nonlinear dynamics due to the mode cou- 
pling between these modes were investigated in Ref. [13] ■ 
Recently, also a coupling of the dipole mode with the 
breathing and quadrupole mode was analyzed in the im- 
mediate vicinity of a Feshbach resonance [ij . 

It is well known that the majority of theoretical stud- 
ies of collective excitations of a BEC is mainly focused 
on considering the two-body contact interaction due to 
the diluteness of the quantum gas [13, [H-il IH, El]. 
On the other hand, in view of the experimental progress 
with BECs in atomic waveguides and on the surface of 
atomic chips, which involve a strong compression of the 
traps, a significant increase of the density of BECs can 
be achieved, thus also three-body interactions may play 
an important ro le 145144 7l |. Theoretical and experimental 
studies m m |49l for a BEC of ^'^Rb atoms indicated 
that the real part of the three-body interaction term can 
be 10"^ — 10* times larger than the imaginary part. The 
imaginary part, which arises from three-body recombi- 
nations, limits the lifetime of the condensate. However, 
even for a small strength of the three-body interaction, 
the region of stability for the condensate can be extended 
considerably according to Ref. [50|. Due to the three- 
body interaction, the density profile of the condensate 
changes [5l| and, correspondingly, also the excitation 
spectrum of the collective oscillations is modified [52l.[53j. 
The modulation instability of a trapped BEC immersed 
in a highly elongated harmonic trap with two- and three- 
body interactions was studied both analytically and nu- 
merically [54]. The effect of the three-body interaction 
was furthermore studied in ultracold bosonic atoms in 
an optical lattice 47, 55 61], BCS-BEC crossover [13], 
complex solitons BEC [63] and vortex BEC [&i ]. 

The collective excitations of a ID BEC in a quadratic 
plus quartic trap are studied on the level of the Gross- 
Pitaevskii equation, and the collective excitation spectra 
are calculated as functions of the two-body interaction 
as well as an anharmonic parameter in Ref. i65i]. The 
collective excitations of a one-dimensional Bose-Einstein 
condensate trapped in an anharmonic potential were also 
studied in Ref. [663 by solving the time-dependent Tonks- 
Girardeau equation. In Ref. [67] the authors investigated 
the collective excitations of a ID BEC with two- and 
three-body interactions in a harmonic trap with an an- 
harmonic distortion up to a quintic term, by mainly fo- 
cusing on the effect of the three-body interaction on the 
excitations by a variational analysis of the (CP) equation. 
Dynamical properties of a rotating BEC confined in an 
anharmonic trap are investigated in Ref. [68] . The transi- 
tion temperature, the depletion of the condensate atoms, 
and the collective excitations of a BEC with two- and 
three-body interactions in an anharmonic trap at finite 
temperature are studied in Ref. . Ref. ^3] shows that 
the frequency of the collective excitation is also signifi- 



cantly affected by the strength of the three-body interac- 
tion and the anharmonicity of the potential. In Ref. [7l| 
the authors investigated the collective excitations and the 
stability of a BEC in a one-dimensional trapping geom- 
etry for the case of a repulsive or attractive three-body 
interaction together with a repulsive two-body interac- 
tions by a standard variational approach. 

An attractive two-body interaction has a profound ef- 
fect on the stability of a BEC, since a large enough attrac- 
tive interaction will cause the BEC to become unstable 
and collapse [i,[72-l73. In Refs. [Hlzi-lzi it was shown 
within the Gross-Pitaevskii formalism, that the critical 
particle number for a BEC in cylindrical traps can be 
obtained numerically. In Ref. [79] the CP equation for a 
system composed of attractive bosons confined in a har- 
monic trap was analyzed via a controlled perturbation 
theory. In this way it was shown that the critical particle 
number strongly depends on the anisotropy of the trap. 
In Ref. JSO] , an analytical solution of the Gross-Pitaevskii 
equation for a BEC with negative effective interaction 
strength was considered in a cylinder-symmetric trap. In 
particular, an analytical formula was derived which an- 
alyzes the critical particle number as a function of the 
anisotropy of the confining potential. Ref. [sTj investi- 
gates how the critical particle number changes due to 
small anharmonic terms added to the confining potential 
of an atomic condensed system with negative two-body 
interaction. The dynamics of a recombination three- 
body Gross-Pitaevskii equation for trapped atomic sys- 
tems with attractive two-body interaction was investi- 
gated numerically [11] • In Refs. [t^,!!^] the authors show 
that, in a dilute g small repulsive three-body interac- 
tion added to an attractive two-body interaction is able 
to stabilize the condensate such that the critical number 
of atoms in traps increases [721 ]. The maximal critical 



number of atoms occurs in a spherical trap according to 
various theoretical predictions [H, [tI, [t^, 113, Isl] , which 
agree with experimental measurements [6|. 

In this paper, we study the dynamics of the condensate 
in general and its collective oscillation modes in partic- 
ular by changing the geometry of the trapping poten- 
tial. The asymmetry of the confining potential leads to 
important nonlinear effects, including resonances in the 
frequencies of collective oscillation modes of the conden- 
sate [22, 2^, [s^] . The BEC dynamics at zero temperature 
is described by the nonlinear GP equation for the con- 
densate wavefunction with two- and three-body interac- 
tions. Within a variational approach the partial differ- 
ential equation of Gross and Pitaevskii is transformed in 
Sec. |TT] into a set of ordinary differential equations for the 
widths of the condensate in an axially-symmetric har- 
monic trap with both two- and three-body interactions. 
Although this represents a significant simplification in 
the description of a BEC, Ref. [23] demonstrates explic- 
itly that these ordinary differential equations yield ap- 
proximatively even for long evolution times the correct 
dynamics of the widths. We discuss in detail in Sec. IIIII 
the resulting stability of the condensate. First, we con- 
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sider the case of an attractive two-body interaction and 
a vanishing three-body interaction. Then we consider 
the case when we have attractive two-body and repul- 
sive three-body interaction. In Sec. IIVI we study in de- 
tail geometric resonances and derive an explicit analytic 
results for the frequency shifts for the case of an axially- 
symmetric condensate with two- and three-body interac- 
tions based on a perturbative expansion and a Poincare- 
Lindstedt method. This frequency shift is calculated for 
a quadrupole mode in Sec. IIV A[ for a breathing mode 
in Sec. IIVBI and the derived analytical results are then 
compared with the results of numerical simulations in 
Sec. IIV CI In that section we also compare results of nu- 
merical simulations for radial and longitudinal widths of 
the condensate and the corresponding excitations spectra 
with the analytical results obtained using perturbation 
theory. Then, in Sec. |Vl we analyze the resonant mode 
coupling and the generation of second harmonics of the 
collective modes, which are induced by nonlinear effects. 
At first we consider a BEC in the initial state correspond- 
ing to the stationary ground state with a small pertur- 
bation proportional to the eigenvector of the quadrupole 
mode, which leads to quadrupole mode oscillations. In 
the linear case, we have small-amplitude oscillations of 
the condensate size around the equilibrium widths, and 
we are in the regime of linear stability analysis. However, 
when the frequencies of collective modes are approached, 
we obtain a resonant behavior which is characterized by 
large amplitude oscillations. In this case it is clear that a 
linear response analysis does no longer provide a qualita- 
tively good description of the system dynamics [23|. Fi- 
nally, in Sec. IVII we summarize our findings and present 
our conclusions. 



tions, respectively. The two-body nonlinearity 172 is pro- 
portional to the s-wave scattering length a, and is given 
by (72 = ^Trfi^a/m, where m denotes the mass of the cor- 
responding atomic species. 

The three-body nonlinearity becomes important for 
large values of the scattering length, but also for small 
values of a close to the ideal gas regime. It is well known 
that the stability against the collapse of *^Rb cannot be 
described by using only the two-body scattering. The 
three-body scattering also plays an essential role in un- 
derstanding the Efimov physics (87|, 188,] • Braaten and Ni- 
eto d^l have used an effective field theory to calculate the 
strength of the three-body interaction which effectively 
arises from the two-body interaction, and obtained the 
resuh 53 (k) = 3847r(47r - 3%/3)[ln Ka + B]h'^a'^ /m, where 
K is an arbitrary wave number and i? is a complex con- 
stant, both being suitably fixed in Ref. The effective 
three-body coupling strength represents a complex num- 
ber, where its imaginary part describes recombination 
effects. However, its real part is much larger, and the 
fit to experimental data for ^^Rb and ^^Rb gives typi- 
cal values for Ke(a£)/h of the order of lO^^^cm^s"^ to 
10-26cm6s-i dilTi US- 
Equation ([1} can be cast into a variational problem, 
which corresponds to the extremization of the action de- 
fined by the Lagrangian 

L{t)^ J C{r,t)dr, (2) 

with the Lagranian density 



II. VARIATIONAL APPROACH 

The dynamics of a Bose-Einstein-condensed gas in a 
trap at zero temperature is well described by the time- 
dependent GP equation [t^, Izll ■ Usually, only two-body 
contact interactions are considered due to the diluteness 
of the gas. In this paper, however, we study systems 
where also three-body contact interactions have to be 
taken into account [48|,|86i]. In that case, the GP equation 
has the form 



— |VV|'-F(r)|V|' 
(3) 



In order to analytically study the dynamics of BEC sys- 
tems with two- and three-body interactions, we use the 
Gaussian variational ansatz, which was introduced in 
Refs. [H, [11]. For an axially symmetric trap, this time- 
dependent ansatz reads 
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A + V{r)+g2N\^{r,t)f 



(1) 



where ^(r, denotes a condensate wave function nor- 
malized to unity, and N is the total number of atoms 
in the condensate. On the right-hand side of the 
above equation we have a kinetic energy term, an exter- 
nal axially-symmetric harmonic trap potential V(r) = 
■^mujp [p^ -I- A^z^) with the anisotropy parameter A = 
ujz/^p, while the parameters 172 and 53 account for the 
strength of two-body and three-body contact interac- 



i/i^(p,z,i) = A^(t)exp 
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where Af{t) = 1/ 
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{t)uz{t) is a normalization fac- 



tor, while Up{t), Uz{t), (t>zit), and 0p(i) are variational 
parameters, representing radial and axial condensate 
widths and the corresponding phases. If we insert the 
above Gaussian ansatz into the Lagrangian (jS]), we ob- 
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tain the Lagrange function 
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From the corresponding Euler-Lagrange equations we ob- 
tain the equations of motion for ah variational parame- 
and (f)z can be expressed exphcitly 

ac- 



ters. The phases 
in terms of 
cording to 



in terms of first derivatives of the widths Up and u 
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Inserting Eqs. ^ into the equations for the widths, we 
get second-order differential equation for Up and Uq. Af- 
ter introducing the dimensionless parameters 

Wi^Wi/wp, Ui^Ui/t, i^t/ujp (7) 



with the oscillator length ^ — ^h/mujp, we obtain a sys- 
tem of two second-order differential equations for Up and 
Uz in the dimensionless form 
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(8) 
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where, for simplicity, we drop the tilde sign in the dimen- 
sionless widths. In the above equations, 

p ^ g2N/{27rf^^hujp£^ = ^/2/^Na/e , (10) 
k = 2g3N^/9V3TT^ujpM^ , (11) 

where p denotes the dimensionless two-body interaction 
strength, while the parameter k is the dimensionless 
three-body interaction strength, which can be also ex- 
pressed in terms of p as 



IQgstK^p 2 



(12) 



For iV = lO'^ atoms of ^'^Rh [40, |43] in a trap with 
ujp = 27r X 112 Hz, the two-body interaction strength 
is g2 ^ 5h X lO^^^cm^s"^, yielding p = 426. The three- 
body interaction is of the order g^ H x 10~26cni^s~^ 
[47j , which gives the dimensionless three-body interaction 
value k = 525. 

Although the value of k is larger than that of p, the 
corresponding terms in Eqs. (|5]) and i.e. k/u^ul and 
k/UpU^, are suppressed by the factor u^Uz compared to 
the respective p-terms. This factor can be estimated by 



taking into account the equilibrium positions Upo and Uzo, 
which are obtained by solving the stationary equations 
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(13) 
(14) 



For the anisotropy A = 3/2, one numerically obtains 
UpQ ~ 3.67 and Uzo ~ 2.45, yielding the value u'^qUzo ~ 
33. This shows that the terms proportional to k have 
the effective coupling fc/33 « 15.9, which makes them 
small corrections of the order of 4% to the leading two- 
body interaction terms. However, if the system exhibits 
resonances, this may no longer be true anymore, and 
three-body interactions can play a significant role for the 
system dynamics. In this paper we study geometric res- 
onances, where it turns out to be necessary to take into 
account effects of three-body interactions. 

Equation ([T]) can be rescaled to a dimensionless form 
with respect to the s-wave scattering length [U Ull, 
m, [t^, [Zli, which can be tuned to any value, large or 
small, positive or negative by applying an external rnag- 
netic field, using the Feshbach resonance technique ^l[. 
Therefore, in this paper we will consider a range of ex- 
perimentally realistic values for dimensionless interaction 
strengths p and k. 

Using the Gaussian approximation enables us to ana- 
lytically estimate frequencies of the low-lying collective 
modes [l^, [H, HI]. This is done by linearizing Eqs. dU 
and © around the equilibrium positions. If we ex- 
pand the condensate widths as Up{t) = UpQ + 5up{t) and 
Uz{t) = Uzo + 5uz{t), insert these expressions into the 
corresponding equations, and expand them around the 
equilibrium widths by keeping only linear terms, we im- 
mediately get the frequencies of the breathing and the 
quadrupole mode 
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and the corresponding right-hand eigenvectors are given 
by: 
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The quadrupole mode has a lower frequency and is 
characterized by out-of phase radial and axial oscil- 
lations, while in-phase oscillations correspond to the 
breathing mode. Figure [T] shows the frequencies of both 



0.0 0.5 1.0 1.5 2.0 2.5 3.0 

A 

FIG. 1: Frequencies of the breathing and quadrupole modes 
versus trap aspect ratio A for p = 1, k — 0.001 (solid lines) 
and p — 10, k — 0.1 (dashed lines). 

modes as functions of the trap aspect ratio A. We see 
that the collective mode frequencies strongly depend on 
the trap anisotropy, whereas a variation of the dimension- 
less interaction strengths p and k yields only marginal 
changes. 



III. STABILITY DIAGRAM 

In this section we discuss the stability of a Bose- 
Einstein condensate in the mean-field framework for sys- 
tems with two- and three-body interactions in an axially- 
symmetric harmonic trap. To study three-body effects, 
we consider several cases of interest: repulsive and attrac- 
tive pure two-body interactions, attractive two-body and 
repulsive three-body interactions, and attractive two- 
and three-body interactions. 

If this system of equations does (not) have bounded 
solutions in the vicinity of the equilibrium widths deter- 
mined by Eqs. ^ and ([9]), then the condensate is consid- 
ered stable (unstable). For the case of a pure repulsive 
two-body interactions, we can immediately see that the 
condensate is always stable. For the case of an attractive 
two-body interaction, the situation is quite different: the 
above system of equations can have no equilibrium, or 
it could have up to three equilibrium solutions. Their 
stability has to be studied by considering small oscilla- 
tions around them. The results of a detailed numerical 
analysis are summarized in Fig. [2] 

The solid line in Fig. [5] represents the critical stability 
line as a function of the trap aspect ratio A for fc = 0. 
Above the solid line, the system has one stable and one 
unstable solution, while below this line the condensate is 
unstable. For A = 0, which corresponds to the limit of 
a cigar-shaped condensate, we have the critical value of 
two-body interactions pc = —0.6204, precisely coinciding 
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FIG. 2: Stability diagram of a BEG as a function of a trap 
aspect ratio A and an attractive dimensionless two-body in- 
teraction strength p. The solid blue curve represents a critical 
stability line for fe = 0, below which the condensate is unsta- 
ble and above which one stable and one unstable equilibrium 
exist. For the case of a small repulsive dimensionless three- 
body interaction strength k = 0.005, the critical stability line 
is given by a doted red line, above which one stable and two 
unstable equilibria exist. In both cases, the dashed black line 
delineates two stable regions: the one with one stable and one 
(for fc = 0) or two (for fc > 0) unstable solutions (designated 
lS-(-unstable), and the one with only one stable solution (des- 
ignated IS). As we can see, the stability region is significantly 
extended for non- vanishing k. 

with Ref. [IB]. For the isotropic case, when A = 1, the 
critical value is Pc ~ —0.535, which again coincides with 
the value from the literature d, [H, [tI, It^] ■ 

Now, if we consider the case of an attractive two-body 
interaction and a small repulsive three-body interaction, 
the results of the stability analysis are qualitatively sim- 
ilar. The system can either have one or three solu- 
tions, and the border between the corresponding regions 
is given by the dashed line. Above the dashed line, only 
one stable real positive solution exists, and the conden- 
sate is stable. Between the dashed and dotted line there 
are three solutions, one being stable and two unstable, 
making the condensate still stable. Quantitatively, we 
see that the critical stability line lies lower, as shown in 
Fig. [21 where the dotted line now presents the critical 
stability line for k = 0.005. 

To further illustrate the above stability analysis, we 
plot in Fig. [3] the frequencies of the low-lying collective 
excitation modes as functions of an attractive two-body 
interaction for the trap anisotropy A — 117/163 The 
solid curves correspond to the case when three-body in- 
teractions are neglected, i.e. k = 0, and we can see that 
the condensate collapses for pc = —0.561. The dashed 
lines give the frequencies for the case when the small 
attractive dimensionless three-body interaction strength 
k — —0.005 exists, and in that case the condensate col- 
lapses already for pc = —0.54. On the other hand, if 
the small repulsive dimensionless three-body interaction 
strength k = 0.005 is present, we see from the dotted lines 
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FIG. 3: Frequencies of the collective excitation modes: 
breathing (B), radial quadrupole (RQ), and quadrupole (Q), 
as functions of an attractive dimensionless two-body interac- 
tion strength p for the trap anisotropy A = 117/163. Solid 
lines correspond to the pure two-body interaction case (k — 
0), dashed (dotted) lines correspond to a small repulsive (at- 
tractive) dimensionless three-body interaction strength k = 
0.005 (fc = -0.005). 



that the stability region of the condensate is extended, 
and the critical value for the collapse is Pc = —0.579. 



IV. FREQUENCIES OF COLLECTIVE MODES 

Close to resonances, the nonlinear structure of the GP 
equation leads to shifts in the frequencies of collective 
oscillation modes compared to the respective values in 
Eq. ([T5|). which are calculated using the linear stability 
analysis. The resonances can be generated by two differ- 
ent mechanisms: they could be induced by an external, 
parametric driving of the system [23| . or emerge due to 
the geometry of the system [ll| . The latter case of geo- 
metric resonances is further investigated in this section. 
In Sec. |V] we show then that geometric resonances also 
induce a coupling of the collective modes, which is only 
possible due to the nonlinear nature of the system dy- 
namics which is due to the interaction. 

In order to analytically study both of these non- 
linearity-induced effects, we apply the standard Poincare- 
Lindstedt method 23., 92-95] in order to develop a per- 
turbation theory that describes dynamics of a BEC with 
two- and three-body interactions. 



A. Quadrupole mode 

We start with working out a perturbation theory for 
the BEC dynamic, which is based on a set of ordinary dif- 
ferential equations ([S])-®, by expanding the condensate 



widths in the series 

Up{t) = Upo + eupi{t) + s'^Up2{t) + s^Up3{t) + ... , (18) 
uz{t) = Uzo + eu,i{t)+e'^u,2{t)+e^u,3{t) + ... ,{19) 

where the smallness parameter e stems from the respec- 
tive initial conditions. Here we study the system dynam- 
ics with the initial conditions in the form 



u(0) = uo + euQ , 
u(0) = 0, 



(20) 
(21) 



when the system is close to the equilibrium position Uq, 
and is perturbed in the direction of the quadrupole os- 
cillation mode eigenvector uq determined by Eq. (jl7l) . 
By inserting expansions (IT^ and (ITOl) into Eqs. ([5])-([n]), 
we obtain the following system of linear differential equa- 
tions: 



{t) + miUpnit) + ■m2Uzri{t) = Xpnit) , (22) 
Uzn{t) +'2m2Upn{t) +m3Uznit) = Xzn{t) , (23) 

where the index n takes integer values n — 1,2,3,.., and 
the quantities mi, rn2, and 7713 are already defined by 
expressions (1161) . The functions Xpn (t) and Xzn (t) depend 
only on the solutions Upi{t) and Uzi{t) of lower orders i, 
i.e. those corresponding io i < n. At first, for n = 1, we 
have xpi(0 — Xzi{t) — and solve the above system of 

equations for Ui(t) — ( "p^^^I 

use the previous solution, and solve the system for U2 (t) , 
where the functions Xp2{t) and Xz2it) are given by 



Then, for n = 2, we 
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Xp2 = aupiu^i + I3u^^ + -fpUp^ , 
Xz2 = 4^UpiM2i -I- +7z'"^i -I- a-Upi 
with the abbreviations 

3p 



13'- 
6 

6 



7p 



Iz 



4 2 






'"pO'^zO 




3k 


3 3 

^pO^zO 


5 4 

^pO'^zO 


6» 


15A: 




7 2 


UpQUzO 


■"pO^zO 


3p 


6k 


2 4 

"p0^20 


"pO^zO 



(24) 
(25) 

(26) 
(27) 
(28) 
(29) 



At each level n of this procedure we use the initial con- 
ditions from Eqs. (UHl) and (PT|) . 

In order to decouple the system of equations 
we use the linear transformation 



Upn{t) = Xn{t) + Vnit) , 
Uzn{t) = CiXn{t) + C2y„(t) 



with the coefficients 
TO3 



Cl,2 



mi =F (m3 — mi)2 + 8m| 



2m2 



(30) 
(31) 



(32) 
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It decouples the system at level n and leads to two inde- 
pendent linear second-order differential equations: 



Xn{t) + UjQXn{t) 



2 u\ , ^■^Xpnit) - Xzn{t) 
'Q"" 

yn{t) +UJ%yr,{t) 



Cl - C2 
Xzn{t) - CiXpn{t) 



Cl ~ C2 



0, (33) 
= . (34) 



From this we conclude that Xn{t) and yn{t) correspond 
to quadrupole and breathing mode oscillations, respec- 
tively. Although the system is perturbed only in the 
direction of the quadrupole mode eigenvector, due to the 
nonlinearity of the system, the breathing mode is excited 
as well. 

The solutions of the above equations depend essen- 
tially on the nature of the inhomogeneous terms, which 
are given by polynomials of harmonic functions of wgi, 
ujBt and their linear combinations {kujQ + mu)B)t- If the 
inhomogeneous part of the equation for Xn {yn) does not 
contain a linear term proportional to the harmonic func- 
tion in uqt {ajst), the particular as well as the homo- 
geneous solutions will be again polynomials in harmonic 
functions of wgt, ajst and their linear combinations. In 
that case, the system exhibits usual small oscillations, 
which are characterized by the frequencies ujq and ujb 
from the linear stability analysis. Compared to linear 
systems, the exception here is that higher harmonics and 
linear combinations of the modes emerge due to the in- 
herent nonlinearity. 

However, if the inhomogeneous part of the equation 
for Xn (yn) contains a linear term proportional to the 
harmonic function in uqI [uj^t), then the corresponding 
particular solution will contain a secular term, i.e. one 
proportional to the t time a harmonic function of ugt 
(uBt), which turns out to be responsible for the frequency 
shift of the collective mode [HI . This happens for the first 
time at level n = 3, when Eqs. can be written 

in vector form as 

U3 (t) + Mu3 (i) -I- Iq,3 cos wgi + . . . = , (35) 

with the matrix 



M : 



mi 7712 

2m2 TO3 



(36) 



and the dots represent the inhomogeneous part of the 
equation, which does not contain linear terms propor- 
tional to harmonic functions in iogt. The expression for 
Ig.n can be calculated systematically for any given level 
n in Mathematica software package [96|. The particular 
solution of Eq. (1351) has the form 



U3.p(i) = 



— £ 



2("q) 



'-Q-3 



2ujQ 



with 



Ugt sin (jjgt - 



2m2 



'4mn 



mi 



mi 



(37) 



(38) 



where Ug denotes the left-hand eigenvector correspond- 
ing to the quadrupole mode and, again, dots represent 
part of the particular solution due to other inhomoge- 
neous terms. The complete solution of Eq. ([55)1 is given 
by the sum of the homogeneous solution uq cos luqI and 
the particular solution u^ pit). The secular term can 
be now absorbed by a shift in the quadrupole mode fre- 
quency according to 



U3{t) = Uq COS UJQt - e 



(u^)^^1q,3 



2ujQ 

Uq COs(wq + ALLlQ)t + 



UQtsinujQt + 



which is quadratic in the smallness parameter e: 



ujq{£) = + AujQ =ujQ-e' 



2ujQ 



(39) 



(40) 



The expression for (uq)-^Iq_3 is most easily calculated in 
Mathematica software package (96j , and reads 



Q 



{lob - '2ujq){ujb + 2a;Q) 



(41) 



where /q^3 is a regular function, without poles for real 
values of its arguments. Therefore, the frequency shift of 
a quadrupole mode to lowest order in e is given by 



2 /q,3 i^Q ,UJB,UpO,UzO,p,k,X) 

2ujq{llib - 2u!q){ujb + Zwq) 



(42) 



and depends explicitly on the trap aspect ratio A, as well 
as implicitly through the A-dependence of the frequen- 
cies loq and ujb- The expression (|42|) for the frequency 
shift has a pole for ujb = 2a;Q. Taking into account the 
fact that ujQ < ujb, as we can see from Eq. (|15p and 
Fig. [l] this condition can, in principle, be satisfied. This 
is designated to be a geometric resonance, since it can be 
obtained just by tuning the geometry of the experiment, 
i.e. the value of the trap aspect ratio A. Higher-order cor- 
rections to AwQ in £ can also be obtained systematically 
by using the developed perturbation theory. In Sec lIVCl 
we compare the analytical result for the frequency shift 
of a quadrupole mode ((42|) with corresponding results of 
numerical simulations. 



B. Breathing mode 

In a similar manner, we study the dynamics of a 
cylindrically-symmetric BEC system when initially only 
the breathing mode is excited. 



u(0) = Uo + £Ub , 
u(0) = 0. 



(43) 
(44) 



Applying again the Poincare-Lindstedt perturbation the- 
ory, we calculate the breathing mode frequency shift. 



UJB (e) = W_B + Aws = ws - £ 



2(u|riB£ 

2ujB 



(45) 
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FIG. 4; A comparison of analytic (solid lines) and numeric 
(dots) results for a nonlinear BEG dynamics with a pure re- 
pulsive two-body interaction p = 1, A: = 0, and e = 0.1. Top 
panels show dynamics of (a) longitudinal and (b) radial con- 
densate widths for the trap aspect ratio A = 2.3 as a function 
of the dimensionless time parameter Upt; bottom panels show 
dynamics of (c) longitudinal and (d) radial condensate widths 
for A = 0.55. 



FIG. 5; A comparison of analytic (solid lines) and numeric 
(dots) results for nonlinear BEG dynamics for a repulsive two- 
body interaction p = 1 and a repulsive three-body interaction 
k — 0.001, with e — 0.1. Top panels show dynamics of (a) lon- 
gitudinal and (b) radial condensate widths for the trap aspect 
ratio A = 0.7 as a function of the dimensionless time parame- 
ter (jjpt; bottom panels show dynamics of (c) longitudinal and 
(d) radial condensate widths for A = 2.05. 



where again the expression (u^)'^Ib.3 is calculated in 
Mathematica software package [96[ . In this way we finally 
yield the following analytic formula for the frequency 
shift of the breathing mode 



2ws(2a;s - a;Q)(2ws -I- wq) 



(46) 



where the function /s,3 is a regular function of its ar- 
guments. Naively looking at this expression, one would 
say that it exhibits a pole for 2ujb — ^q- However, from 
Eq. (IT5|) and Fig. [T]we see that ujq < ujb, and, therefore, 
the condition 2ujb = loq is never satisfied. In the next 
section we numerically demonstrate , that a geometric 
resonance does not occur, and verify the analytical re- 
sult for the frequency shift of the breathing mode. 



parameter Upt for different values of the trap aspect ratio 
A. Analytical solutions are calculated using the third- 
order perturbation theory developed in Sec. IIV Al We 
can see in both figures that the agreement is excellent, 
not only for non-resonant values of the trap aspect ra- 
tio A (top panels), but also for values close to resonances 
(bottom panels), which are characterized by a violent dy- 
namics. Such behavior is just a tell-tale of a geometric 
resonance, and a subsequent analysis of frequency shifts 
is the only proper way to identify these resonances. 

However, before we present this analysis, we show in 
Fig.[6]excitation spectra of BEC the dynamics which cor- 
responds to the initial conditions (P(I1) - (PT|) for p = 1, 
k — 0.001, and two values of the trap aspect ratio, 
A = 1.9 and A = 0.5. For the parameter values in 
Fig. [6ja), the linear stability analysis yields breathing 



C. Comparison with numerical results 

In order to verify our analytical results, we perform 
extensive numerical simulations [97- 100] . At first we fo- 
cus on a description of the nonlinear BEC dynamics, and 
compare our analytical results for the radial and longitu- 
dinal widths of the condensate obtained perturbatively 
to the direct numerical solutions of Eqs. (l8])-([9]). To 
this end we consider a BEC in the initial state corre- 
sponding to the perturbed equilibrium position, where 
the small perturbation is proportional to the eigenvector 
of the quadrupole mode, Eqs. (I^ -(PT |) . Examples of 
the condensate dynamics are shown in Fig. 2] for a pure 
two-body interaction p — \, k = Q with e = 0.1, and in 
Fig.[S]forp = 1, fc = 0.001, e = 0.1. 

In both figures we plot analytical and numerical solu- 
tions for Uz and Up as functions of the dimensionless time 
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FIG. 6: Fourier spectra of the nonlinear BEC dynamics for a 
repulsive two-body interaction p = 1, a repulsive three-body 
interaction k = 0.001, and e = 0.1 for (a) A = 1.9 and (b) A = 
0.5. Each graph shows spectra of both longitudinal and radial 
condensate widths. The locations of all peaks are identified 
as linear combinations of the quadrupole and the breathing 
mode frequency, in correspondence with the analysis based on 
the developed perturbation theory. 
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and quadrupole mode frequencies Eq. (jl5|) with ujb = 
3.65018 and ujq — 1.96208, while for the parameters in 
Fig. Mh) we obtain ub = 2.0084 and ujq = 0.907393 
all expressed in units of Up. In both graphs we can see 
that the Fourier spectra contain two basic modes, loq and 
cjb, whose values agree well with those obtained from the 
linear stability analysis in Eq. (jlSp . and a multitude of 
higher-order harmonics, which are linear combinations of 
the two modes, as pointed out in Sec. IIV Al 

Now we compare the derived analytical results for the 
frequency shifts of the quadrupole and the breathing 
mode with the results of numerical simulations for the 
BEC systems with two- and three-body contact interac- 
tions in a cylindrical trap. In particular, we focus on the 
nonlinear regime and calculate frequency shifts close to 
geometric resonances, which we identify from the poles 
of these shifts. In Figs. [7] and [S] we present the compar- 
ison of analytic (solid lines) and numeric (dots) values 
of relative frequency shifts as functions of the trap as- 
pect ratio A. The analytical results are calculated from 
Eq. ([32]) and respectively, while the numerical data 
are obtained from a Fourier analysis of the excitation 
spectrum, i.e. for each value of A we have calculated the 
corresponding Fourier spectra, as in Fig. [SI and then ex- 
tracted from them the frequency values of the quadrupole 
and the breathing mode. 

In Fig. [TJJa) we show a special case of a pure two-body 
interaction, when fc = 0. The condition for a geomet- 
ric resonance lob = '^.ujq yields the trap aspect ratios 
Al = 0.55 and A2 = 2.056, in good agreement with the 
numerical data, as we can see from the graph. The exis- 
tence of a geometric resonance at Ai = 0.55 is responsible 
for a violent dynamics seen in the bottom panels of Fig.[3J 
as we have pointed out earlier. However, by analyzing the 
frequency shifts we can conclusively show that, indeed, 
the geometric resonance is present. In further graphs we 
see that the excellent agreement between analytical and 
numerical results holds also for other values of p and fc, 
including the case of an attractive two-body interaction 
p — —0.2, which is still within the BEC stability region. 
It is interesting to note the observation that the asymp- 
totic approach to geometric resonances for the case of 
an attractive two-body interaction is reversed compared 
to the case of a repulsive two-body interaction. For in- 
stance, we can see in Fig.[71[d) that Aljq/ujq — > 00 when 
A — ?> A^, and Aujq/ujq — ^ —00 when A — >■ Aj, while 
for an attractive p — —0.2 in Fig. [TKf) we see that the 
situation is reversed. 

In Fig. [5] we compare analytic and numeric results for 
a frequency shift of the breathing mode. As for the 
quadrupole mode, the agreement is excellent for both re- 
pulsive and attractive two-body interaction. As pointed 
out in Sec. IIVBI there are no geometric resonances for 
the breathing mode frequency, since the corresponding 
condition ujq = Iujb cannot be satisfied. 

Finally, we compare our derived analytic results with 
those from Ref. z^, where frequency shifts of collective 
modes were calculated in the Thomas-Fermi (TF) limit 
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FIG. 7: Relative frequency shift of the quadrupole mode as a 
function of the trap aspect ratio A for e = 0.1 and different 
values of two-body and three-body interaction strengths: (a) 
p = 1, = 0, (b) p = 1, fc = 0.001, (c) p = 0.1, fc = 0.001, 
(d) p = 0.1, fc = 0.1, (e) p = -0.2, fc = 0, (f) p = -0.2, 
fc = 0.005. Solid lines represent the analytical result obtained 
from Eq. (|42p . while dots are obtained by a numerical analysis 
of the corresponding excitation spectrum for each value of A, 
as described in Fig. [Sj 
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FIG. 8: Relative frequency shift of the breathing mode as a 
function of the trap aspect ratio A for e = 0.1 and different 
values of two-body and three-body interaction strengths: (a) 
p = 1, fc = 0, (b) p = 0.2, fc = 0.001, (c) p = 0.1, fc = 0.1, 
(d) p — —0.2, fc = 0.005. Solid lines represent the analytical 
result obtained from Eq. (|46p . while dots are obtained by a 
numerical analysis of the corresponding excitation spectrum 
for each value of A, as described in Fig. [6] 



using a hydrodynamic approach. In terms of our vari- 
ational approach, the TF limit corresponds to the limit 
p — ?► 00, so that Eq. (ITSI) for the frequencies of the breath- 
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FIG. 9: Comparison of the analytical results for the relative 
frequency shifts of the quadrupole and the breathing mode in 
the Thomas- Fermi limit from Ref. 22] (solid lines) and the an- 
alytical results derived here using the Poincare-Lindstedt per- 
turbation theory (dashed lines). Panels (a) and (c) present a 
direct comparison of the analytical results derived in Ref. [12] 
using the parabolic variational ansatz and our perturbation 
results derived using the Gaussian ansatz (Q, while panels 
(b) and (d) present the corresponding comparison when the 
parabolic variational ansatz is used in both schemes. 



ing and the quadrupole mode reduces to 



= 2 + ^A^ ± i Vl6 - 16A2 + 9A4 . 



(47) 



in agreement with Ref. [22| . The condition for a ge- 
ometric resonance iob = Swq yields trap aspect ra- 
tios Ai,2 = (7125 ± 729)/ v%, i.e. Ai « 0.683 and 
A2 « 1.952. 

Fig. ini gives a comparison of the relative frequency 
shifts in the TF limit calculated in Ref. [12] using a 
hydrodynamic approach, and our analytical results ob- 
tained using the Poincare-Lindstedt perturbation theory. 
Panels (a) and (c) give a direct comparison, and we see 
that the agreement is good in general. However, we ob- 
serve small differences, which are due to the fact that 
Ref. [l^l uses a parabolic variational ansatz, while we 
use the Gaussian ansatz in Eq. ^ . In panels (b) and (d) 
we compare the relative analytical results when both ap- 
proaches are applied to a system of equations based on a 
parabolic variational ansatz. As expected, the agreement 
is perfect in this case. 



V. RESONANT MODE COUPLING 

In this section we study the phenomenon of non- 
linearity-induced resonant mode coupling. As already 
pointed out, even if a BEG system is excited precisely 
along the quadrupole or, equivalently, the breathing 
mode, the emerging dynamics will lead to small oscilla- 
tions which initially involve only the corresponding mode, 
but then the other collective mode will eventually step in. 



as well as higher harmonics of the two modes and their 
linear combinations. This phenomenon is caused by the 
nonlinear nature of the system [2^. If the trap confine- 
ment of the system allows a geometric resonance, this 
could greatly enhance the mode coupling and speed up 
the emergence of those modes which are initially not ex- 
cited, and therefore we designate it as a resonant mode 
coupling. 

To demonstrate this phenomenon, we use the pertur- 
bative solution of Eqs. ([Hl)-(|ni) with the initial conditions 
defined by Eqs. ([20]) and ([^T|) . for which the initial state 
coincides with the equilibrium with a small perturbation 
proportional to the quadrupole mode eigenvector. The 
second-order perturbative solution can then be written 
in the form 



uo 



ApB 



cos UJst 



(48) 



where dots represent higher harmonics and the respective 



(49) 
(50) 
(51) 
(52) 



Note that the absence of terms linear in e in expressions 
for ApB and A^b is due to the initial condition, i.e. the 
fact that, initially, we only excite the quadrupole mode, 
the constants ci,2 in the above expressions are defined by 
Eq. (|32|). while ApQ2 and ApB2 are calculated to be 

. C27p + ciC2a + c^C2/3 - a - 4ci/3 - c^72 
-^pQ2 = — 7^; , (ooj 



are j 


;iven by 






= £UpQ + 


£ ApQ2 — 2- , 


A^Q 


= ciApQ , 




ApB 


= e'^ApB2 




A^B 


= C2ApB ■ 





A 



pB2 



3(ci - C2) 

-c\P + a - cijp + 4ci/3 - cja + cjjz 

Cl - C2 



(54) 



with a, 13, 7p, and 7^ defined by Eqs. ([^ - (^5)1 . 

In Fig. [To] we see the comparison of the derived analyt- 
ical results and numerical simulations for the amplitudes 
of the breathing mode, which emerge in the second order, 
as expected, due to nonlinearity features of the system. 
The numerical results are obtained, as before, by extract- 
ing the amplitude of the breathing mode from the Fourier 
excitation spectra of the system for each value of the trap 
aspect ratio A. The agreement is quite good, and we see 
again a resonant behavior, which occurs at the same trap 
aspect ratios as for the frequency shift of the quadrupole 
mode. From Eqs. ([51]) and (|52|) we actually see that the 
resonances occur when the condition ujb = 2wq is sat- 
isfied, which is precisely the same condition as for the 
frequency shift. Therefore, geometric resonances are not 
only refiected in the resonances of frequency shifts of col- 
lective modes, but also in the resonant mode coupling. 



11 





(d) 


P 


M 




k 


= 0.001 \ 














0.5 


1 1.5 2 2.5 



(f) 


/7 = -0.2 








k = 0.005 





















I 1.5 
A 



FIG. 10; Amplitudes of the breathing mode emerging in 
the second order of the perturbation theory from nonhnear 
BEG dynamics after initially only the quadrupole mode is 
excited, given as functions of the trap aspect ratio A for 
e — 0.1 and different values of two-body and three-body in- 
teraction strengths p and k. The amplitudes ApB and Azb 
from Eqs. (I51[l and (I52p . corresponds to the longitudinal con- 
densate width of the emerging breathing mode dynamics. 



In addition to the absolute values of the breathing 
mode amplitudes, which are excited through the nonlin- 
ear mode coupling, it is also interesting to look at their 
relative values, compared to the quadrupole mode ampli- 
tudes, i.e. 



R„ — 



A 



pB 



- 2w, 



Q 



2^Q 



A, 



(55) 
(56) 



Figure [TT] shows the comparison of analytical and numer- 
ical results for the relative amplitudes of the nonlinearity- 
excited breathing mode. Due to the geometric reso- 
nances, we see that the trap aspect ratio can be tuned in 
such a way that the resonant mode coupling excites the 
breathing mode with an amplitude which is far larger 
than that of the quadrupole mode, which served as the 
source of excitation. 

Furthermore, from Eqs. (ISTj) and ([5^ we see that, if 
the geometry of the trap is tuned such that ojb = ^^qV^, 
then the amplitudes of the breathing mode vanish simul- 
taneously, i.e. ApB — AzB — 0. Although this is true 
only in the second-order perturbation theory, it still rep- 
resents a significant suppression of the nonlinear mode 
coupling. Therefore, the tunability of the trap aspect 
ratio offers a unique tool for enhancing and suppressing 
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FIG. 11: Ratios of breathing and quadrupole mode ampli- 
tudes emerging in the second order of the perturbation the- 
ory after initially only the quadrupole mode is excited, given 
as functions of the trap aspect ratio A for £ = 0.1 and differ- 
ent values of two-body and three-body interaction strengths p 
and k. The quantity Rp and Rz from Eqs. (|55p and 1)56^ cor- 
responds to the ratio of amplitudes of the breathing and the 
quadrupole mode in the radial and longitudinal condensate 
width. 



the mode coupling in a BEC, which might be of broad 
experimental interest. 

In a similar way, we can initially excite only the breath- 
ing mode, which corresponds to Eqs. ([H])-® with initial 
conditions defined in Eqs. and (H^ . The solution in 
the second-order perturbation theory has again the form 



UO + 1 4''^ ) COSWst + ( ^''^ ) COSWgt 



but now the respective amplitudes read 



ApB 


= eupB + 


c2 A 

£ ■^pB2 — — 1 


A^B 


= C2ApB , 




ApQ 


= £^ApQ2 




AzQ 


= ciApQ , 





(57) 

(58) 
(59) 
(60) 
(61) 



and the coefficients ApB2 and ApQ2 are given by 



- -ci7p - ciC2a - cicl/S + a + 4c2l3 + c^jz 
^pB2 = TTT. 71 ' ^^^i 



A 



pQ2 



3(ci - C2) 

cf /3 - g + C2lp - 4c2^ + cla - cl-f^ 

Cl - C2 



(63) 



12 



0. 12 


(a) 


P 


= 1 




0.1 




k 


= 




0.08 








y' m 


nl" 0.06 










0.04 










0.02 



























0.5 


1 1.5 

A 


2 


0.12 


(c) 


P 


= 1 




0.1 


k 


= 0.001 




0.08 










a;^ 0.06 










0.04 










0.02 






















0.1 
0.08 
0.06 
0.04 
0.02 




0.1 
0.08 
0.06 
°* 0.04 
0.02 




(b) 



p=\ 

k = 



(d) 



0.5 1 

A 

k = 0.001 




FIG. 12: Ratios of quadrupole and breathing mode ampli- 
tudes emerging in the second order of the perturbation the- 
ory after initially only the breathing mode is excited, given 
as functions of the trap aspect ratio A for e = 0.1 and differ- 
ent values of two-body and three-body interaction strengths p 
and k. The quantity Rp and Rz from Eqs. (|64p and 1)65^ cor- 
responds to the ratio of amplitudes of the breathing and the 
quadrupole mode in the radial and longitudinal condensate 
width. 



In this case, the ratios of amphtudes are given by 

,2 



Ro 



ApQ ^ 2^1 

ApB 



4^1 
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zB 



4-1 
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Figure [T^] compares analytical and numerical results 
for the mode coupling when initially only the breathing 
mode is excited, and then the quadrupole mode emerges 
due to the nonlinearity of the system. As for the case 
of the breathing mode frequency shift, there are no res- 
onances, since ljb > i^q, and the resonance condition 
2ujb = wq cannot be satisfied, as we see from the graphs. 
Also, the condition ojbV^ — ujq cannot be satisfied, and 
therefore the amplitude of the quadrupole mode cannot 
be fully suppressed. For a repulsive two-body interaction 
in Fig.[T2][a)-(d) we see that the ratios Rp and Rz are be- 
low 10%, and the mode coupling mechanism is not able 
to produce a significant amplitude for the quadrupole 
mode. For the case of an attractive two-body interaction 
in Fig. [T2l e)-(f). the ratio increases and the generated 
quadrupole mode amplitude is stronger. Here the agree- 
ment between analytical and numerical results is only 
qualitative, so that the perturbation theory has to be 



carried out to higher orders in the small parameter e in 
order to improve the agreement. 



VI. CONCLUSIONS 



We have studied the dynamics and collective excita- 
tions of a BEC for different trap aspect ratios at zero 
temperature. In particular, we have investigated promi- 
nent nonlinear effects that arise due to two- and three- 
body interactions, as well as their delicate interplay. We 
have discussed in detail the stability of a condensate 
in an axially-symmetric harmonic trap for the experi- 
mentally most relevant setups: repulsive and attractive 
two-body interactions, attractive two-body and repulsive 
three-body interactions, and attractive two- and three- 
body interactions. We have shown that even a small 
repulsive three-body interaction is able to extend the sta- 
bility region of the condensate beyond the critical number 
of atoms when two-body interaction is attractive. 

Using a perturbation theory and a Poincare-Lindstedt 
analysis of a Gaussian variational approach, we have 
studied in detail the interplay between nonlinear effects 
due to two- and three-body interactions, and the effects 
of the trap geometry. Within the variational approach 
and the Poincare-Lindstedt method, we have analyti- 
cally calculated frequency shifts and a mode coupling in 
order to identify geometric resonances of collective os- 
cillation modes of an axially-symmetric BEC. We have 
also shown that the observed geometric resonances can 
be suppressed if two- and three-body interactions are ap- 
propriately fine-tuned. 

To verify analytically the derived results, we have 
used extensive numerical simulations, which provide de- 
tailed excitation spectra of the BEC dynamics. We have 
numerically observed and analytically described several 
prominent nonlinear features of BEC systems: signifi- 
cant shifts in the frequencies of collective modes, gen- 
eration of higher harmonics and linear combinations of 
collective modes, as well as resonant and non-resonant 
mode coupling. The inherent nonlinearity of BEC sys- 
tems couples collective modes, and even if the system is 
excited so that the perturbation corresponds initially to 
the eigenvector of a particular mode, the nonlinear dy- 
namics of the condensate will eventually excite also other 
modes, due to the mode coupling. The presence of geo- 
metric resonances can significantly enhance this effect, as 
we have shown using the developed perturbation theory. 
All obtained analytical results are found to be in good 
agreement with the numerical results. Furthermore, the 
results for frequency shifts are shown to coincide with the 
earlier derived analytical results "2^ within the hydrody- 
namic approach in the Thomas-Fermi approximation. 
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